
// Copyright (c) Gong Chengming, NARI-TECH, 2011

#include <iostream>
#include "./lu.h"

int main() {
  typedef double value_t;
  typedef Lu<value_t>::MT matrix_t;

  const int N = 4;

  matrix_t A;
  A.SetDim(N);

  value_t v[N*N] = {30, 1, 3, -5, -2, 12, 9, 4, 8, -2, -70, 3, -1, 3, 22, 22};
  A.Assign(&v[0]);
  std::cout<<"A = \n"<<A<<"\n";

  Lu<value_t> lu(A);

  value_t b1[N] = {0, 1, 0, 0};
  value_t b2[N] = {3, 7, -1, 5};

  lu.Solve(&b1[0]).Solve(&b2[0]);

  Matrix<value_t> rhs1(N, 1);
  rhs1.Assign(b1);
  Matrix<value_t> rhs2(N, 1);
  rhs2.Assign(b2);
  std::cout<<"Result 1: \n"<<rhs1<<"\n Result 2:\n"<<rhs2<<"\n";

  return 0;
}
